Krill Life Parammetrs estimation in a context spatial in WAP
Supporting analysis to incorporate in Krill stock assessment frameworks
Mardones, M; Cárdenas, C.
16 May, 2024
rm(list = ls())
knitr::opts_chunk$set(echo = TRUE,
message = FALSE,
warning = FALSE,
fig.align = 'center',
fig.width = 8,
fig.height = 6,
dev = 'jpeg',
dpi = 300)
#XQuartz is a mess, put this in your onload to default to cairo instead
options(bitmapType = "cairo")
# (https://github.com/tidyverse/ggplot2/issues/2655)
# Lo mapas se hacen mas rapidolibrary(here)
library(kableExtra)
library(ggthemes)
library(ggrepel)
library(ggridges)
#analisis
library(ggpubr)
library(easystats) # multiples unciones analiticas
library(sf)
library(tidyverse, quietly = TRUE)
library(modelsummary)
library(terra) # replace raster
library(TropFishR)
library(mixR)
library(CCAMLRGIS)Background
Methodology
A important piece of information for a stock evaluation refers to the biological components such as average sizes and weights across areas and years. To do this, we will explore the biological data and prepare the output to add it into stock assessment integrate model (Methot & Wetzel, 2013).
The object ohbio2 come from data exploration analysis in
data request CCAMLR data. This objetc have bio information from
krill.
#cargo objeto
meta <- get("METADATA")
c1 <- get("C1")
ohbio <- get("OBS_HAUL_BIOLOGY")Join data set with master as c1 set. This join is
trought obs_haul_id variable to get geoposition
variables
ohbio2 <- left_join(c1, ohbio, by="obs_haul_id")
names(ohbio2)## [1] "c1_id.x" "obs_haul_id"
## [3] "obs_logbook_id.x" "obs_haul_number"
## [5] "haul_number.x" "vessel_name"
## [7] "vessel_nationality_code" "fishing_purpose_code"
## [9] "season_ccamlr" "target_species"
## [11] "asd_code" "trawl_technique"
## [13] "catchperiod_code" "date_catchperiod_start"
## [15] "datetime_set_start" "datetime_set_end"
## [17] "datetime_haul_start" "datetime_haul_end"
## [19] "datetime_timezone" "depth_gear_set_end_m"
## [21] "depth_gear_haul_start_m" "depth_bottom_set_end_m"
## [23] "depth_bottom_haul_start_m" "latitude_set_end"
## [25] "longitude_set_end" "latitude_haul_start"
## [27] "longitude_haul_start" "gear_type_code"
## [29] "gear_type" "mesh_code"
## [31] "trawl_net_number" "notes"
## [33] "trawl_duration_depth_h" "trawl_duration_total_h"
## [35] "krill_greenweight_kg" "c1_id.y"
## [37] "obs_logbook_id.y" "haul_number.y"
## [39] "taxon_code" "taxon_scientific_name"
## [41] "taxon_family" "maturity_stage"
## [43] "sex_code" "length_total_cm"
## [45] "greenweight_kg"Firsts glance. Test how many register have by year. In this case,
length_total_cm by season ccamlr. Same exercise in date
period date_catchperiod_start to separate dates.
ohbio3 <- ohbio2 %>%
mutate(Year = year(date_catchperiod_start),
Month = month(date_catchperiod_start),
Day = day(date_catchperiod_start)) %>%
#toupper() para convertir los valores a mayúsculas
mutate(sex_code = toupper(sex_code))Save data further analysis
length481 <-ohbio3 %>%
dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47, 43, 45) |> filter(asd_code=="481")
#save(length481, file = "length481.RData")First thing is get different rater layer to join krill data length according different porpoises.
# Cargo linea de costa
coast <- load_Coastline()
coast1<- st_as_sf(coast)
coast2 = st_transform(coast1, "+proj=latlong +ellps=WGS84")
# Uso las agrupaciones de Strata
strata <- st_read("~/DOCAS/Mapas/Antarctic_SHPfiles/Strata.shp",
quiet=T)
strata=st_transform(strata, "+proj=latlong +ellps=WGS84")Test Strata SSMU, just to know another way to join data, but this kind of spatial structuration is deprecated to mamagemente use (Figure@ref(fig:maptest2).
Grouping bio data into stratas
names(length481)## [1] "vessel_nationality_code" "season_ccamlr"
## [3] "asd_code" "trawl_technique"
## [5] "date_catchperiod_start" "latitude_set_end"
## [7] "longitude_set_end" "gear_type"
## [9] "maturity_stage" "length_total_cm"
## [11] "Year" "Month"
## [13] "sex_code" "greenweight_kg"ohbio6 <- st_as_sf(length481 %>%
drop_na(latitude_set_end),
coords = c("longitude_set_end", "latitude_set_end"),
crs = "+proj=latlong +ellps=WGS84")Study Area
# y testeo el mapa
ssmap <- ggplot()+
geom_sf(data = strata |>
filter(ID != "Outer") |>
mutate(ID = str_replace_all(ID, "Extra", "GS")),color="red")+
geom_sf(data = ohbio6 |>
drop_na(length_total_cm) |>
filter(Year>2017),
aes(color = length_total_cm)) +
geom_text_repel(data = strata |>
filter(ID != "Outer") |>
mutate(ID = str_replace_all(ID, "Extra", "GS")),
aes(x = Labx, y = Laby,
label = ID),
min.segment.length = 0,
box.padding = 2,
max.overlaps = 10)+
geom_sf(data = coast2, colour="black", fill=NA)+
scale_color_viridis_c(option = "G",
name="Length (cm)")+
ylim(230000, 2220000)+
xlim(-3095349 , -1858911)+
# coord_sf(crs = 32610)+ #sistema de prpyecccion para campos completos
coord_sf(crs = 6932)+
theme_bw()
ssmap
Strata Maps in 48.1
Length composition by Strata CCAMLR to visualization first. First step is group data into to poligons strata.
strata <- st_make_valid(strata)
sf4 <- st_join(strata, ohbio6)
# Save an object to a file
saveRDS(sf4, file = "sf4.rds")Load RData
# Restore the object
sf4 <- readRDS("sf4.rds")Data
histogram length data to viz in anthor way.
jzstrata <- ggplot(sf4 %>%
mutate(ID = if_else(ID == "Extra", "GERLASHE", ID)) %>%
filter(Year>2010,
ID !="Outer"),
aes(x=length_total_cm,
y = as.factor(Month),
fill=ID))+
geom_density_ridges(stat = "binline", bins = 30,
scale = 1.9,
draw_baseline = FALSE,
alpha=0.9)+
facet_grid(Year~ID) +
geom_vline(xintercept = 3.6, color = "red")+
scale_x_continuous(breaks = seq(from = 1, to = 10,
by = 2))+
scale_y_discrete(breaks = seq(from = 1,
to = 12, by = 4))+
scale_fill_viridis_d(name="Strata",
option="F")+
theme_few()+
xlab("Length (cm.)")+
ylab("")
jzstrata
Now, must filter the DF
sf5 <- sf4 |>
dplyr::select(c("ID",
"date_catchperiod_start",
"length_total_cm",
"Year",
"Month",
"sex_code",
"greenweight_kg")) |>
mutate(ID = if_else(ID == "Extra", "GERLASHE", ID)) |>
filter(ID !="Outer") |>
data.frame()Calculating the proportion of records per year and per stratum in a
table. This table allows us to establish criteria for defining the
strata on which we will analyze krill parameter. We leave out just
JOIN.
# Calcular la proporción
proporcion <- (round(table(sf5$ID) / sum(table(sf5$ID))*100,2))
# Crear una tabla con las proporciones
tabla_proporcion <- as.data.frame(proporcion)
tabla_proporcion$Categoria <- rownames(tabla_proporcion)
rownames(tabla_proporcion) <- NULL
# Renombrar las columnas
colnames(tabla_proporcion) <- c("Strata", "%")
# Mostrar la tabla utilizando kbl()
kbl(tabla_proporcion,
caption = "Proportion length register by Strata") |>
kable_classic(full_width = F,
html_font = "Cambria") |>
kable_styling(bootstrap_options = "striped",
latex_options = "striped")| Strata | % | NA |
|---|---|---|
| BS | 62.71 | 1 |
| EI | 6.58 | 2 |
| GERLASHE | 12.13 | 3 |
| JOIN | 0.14 | 4 |
| SSIW | 18.43 | 5 |
Parameters estimation k and L_inf_
Method ELEFAN ((Tobias2017?))
sf5fil <- sf5 |>
filter(ID %in% c("BS", "EI", "GERLASHE", "SSIW")) |>
mutate(date_catchperiod_start = as.Date(date_catchperiod_start)) |>
mutate(yearly_Group = floor_date(date_catchperiod_start, "year")) |> # New column
drop_na(length_total_cm)
# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLASHE", "EI")
# Crear una lista para almacenar los resultados de cada ID
lfq_results <- list()
# Iterar sobre cada ID
for (id in unique(sf5fil$ID)) {
# Obtener el nombre correspondiente
nombre <- nombres[which(unique(sf5fil$ID) == id)]
# Filtrar los datos para el ID actual
df_id <- sf5fil %>% filter(ID == id)
# Crear un objeto lfq para el ID actual
lfq <- lfqCreate(data = df_id,
Lname = "length_total_cm",
Dname = "yearly_Group",
bin_size = 0.1)
# Agregar el resultado a la lista
lfq_results[[id]] <- lfq
# Graficar el objeto lfq
plot(lfq, Fname = "catch",
main = nombre)
}
Identified distribution
ei <- Bhattacharya(lfq_results$EI)## Interactive session needed for Bhattacharya.ss <- Bhattacharya(lfq_results$SSIW)## Interactive session needed for Bhattacharya.bs <- Bhattacharya(lfq_results$BS)## Interactive session needed for Bhattacharya.gs <- Bhattacharya(lfq_results$GERLASHE)## Interactive session needed for Bhattacharya.ahora damos asignaciones a curvas para representación
# Crear una lista para almacenar los resultados de PW_results <- list()
# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLASHE", "EI")
# Iterar sobre cada objeto lfq almacenado en lfq_results
for (i in seq_along(lfq_results)) {
lfq <- lfq_results[[i]]
# Obtener el nombre correspondiente
nombre <- nombres[i]
# Restructurar el objeto lfq
lfqbin <- lfqRestructure(lfq, MA = 3, addl.sqrt = TRUE)
# Graficar el objeto lfq reestructurado
plot(lfqbin, hist.col = c("white", "black"),
image.col = c(rep(rgb(1,0.8,0.8),1000), "white",
rep(rgb(0.8,0.8,1),1000)),
ylim = c(0,max(lfqbin$midLengths+0.5)),
main = nombre)
# Ajustar curvas al objeto lfq reestructurado
tmp <- lfqFitCurves(lfqbin, par = list(Linf=6.5,
K=0.45,
t_anchor=0.5),
draw = TRUE, col=4, lty=2)
}
Proceed to fit Modal Progression analysis with
ELEFAN
Parameters to SSIW
# Response surface analyss
res_RSAsswi <- ELEFAN(lfq = lfq_results$SSIW,
MA = 5,
Linf_range = seq(4.5, 7,
length.out = 30),
K_range = exp(seq(log(0.1),
log(2),
length.out = 30)),
C = 0.5,
ts = 0.5,
method = "cross",
cross.date = lfq_results$SSIW$dates[3],
cross.midLength = lfq_results$SSIW$midLengths[5],
contour = TRUE,
add.values = FALSE,
agemax = 5,
hide.progressbar = TRUE)## Optimisation procuedure of ELEFAN is running.
## This will take some time.
## The process bar will inform you about the process of the calculations.# run ELEFAN with simulated annealing
res_SAsswi <- ELEFAN_SA(lfq_results$SSIW,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 4.5,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 4,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.011048 maxTime = 30
## Emini is: -0.2661270789
## xmini are:
## 5.328667012 0.8705196602 0.2853670418 0.4288634807 0.1844906621
## Totally it used 30.011068 secs
## No. of function call is: 1092# run ELEFAN with genetic algorithm
res_GAsswi <- ELEFAN_GA(lfq_results$SSIW,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 4.5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
Table with different method
# show results
RSAsswi <- unlist(res_RSAsswi$par)
GAsswi <- unlist(res_GAsswi$par)
SAsswi <- unlist(res_SAsswi$par)
parsswi <- round(rbind(RSAsswi, GAsswi, SAsswi),3)
parsswi %>%
kbl(booktabs = T,
position="ht!",
caption = "Parametres LH to SSWI") %>%
kable_paper("hover",
full_width = F)%>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE,
font_size=9)| Linf | K | t_anchor | C | ts | phiL | |
|---|---|---|---|---|---|---|
| RSAsswi | 5.448 | 0.640 | 0.790 | 0.500 | 0.500 | 1.279 |
| GAsswi | 5.879 | 0.596 | 0.331 | 0.463 | 0.475 | 1.314 |
| SAsswi | 5.329 | 0.871 | 0.285 | 0.429 | 0.184 | 1.393 |
Parametrs to EI
# Response surface analyss
res_RSAei <- ELEFAN(lfq = lfq_results$EI,
MA = 5,
Linf_range = seq(4.5, 7,
length.out = 30),
K_range = exp(seq(log(0.1),
log(2),
length.out = 30)),
C = 0.5,
ts = 0.5,
method = "cross",
cross.date = lfq_results$EI$dates[3],
cross.midLength = lfq_results$EI$midLengths[5],
contour = TRUE,
add.values = FALSE,
agemax = 5,
hide.progressbar = TRUE)## Optimisation procuedure of ELEFAN is running.
## This will take some time.
## The process bar will inform you about the process of the calculations.# run ELEFAN with simulated annealing
res_SAei <- ELEFAN_SA(lfq_results$EI,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 4.5,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 4,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.035548 maxTime = 30
## Emini is: -0.2111109108
## xmini are:
## 5.863710068 0.9813904597 0.873453021 0.9993372839 0.5819024704
## Totally it used 30.035574 secs
## No. of function call is: 1329# run ELEFAN with genetic algorithm
res_GAei <- ELEFAN_GA(lfq_results$EI,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 4.5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
# show results
RSAei <- unlist(res_RSAei$par)
GAei <- unlist(res_GAei$par)
SAei <- unlist(res_SAei$par)
parei <- round(rbind(RSAei, GAei, SAei),3)parei %>%
kbl(booktabs = T,
position="ht!",
caption = "Parametres LH to Elephand Island") %>%
kable_paper("hover",
full_width = F)%>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE,
font_size=9)| Linf | K | t_anchor | C | ts | phiL | |
|---|---|---|---|---|---|---|
| RSAei | 5.879 | 0.470 | 0.750 | 0.500 | 0.500 | 1.211 |
| GAei | 5.737 | 0.595 | 0.960 | 0.618 | 0.528 | 1.292 |
| SAei | 5.864 | 0.981 | 0.873 | 0.999 | 0.582 | 1.528 |
Parametrs to BS
# Response surface analyss
res_RSAbs <- ELEFAN(lfq = lfq_results$BS,
MA = 5,
Linf_range = seq(4.5, 7,
length.out = 30),
K_range = exp(seq(log(0.1),
log(2),
length.out = 30)),
C = 0.5,
ts = 0.5,
method = "cross",
cross.date = lfq_results$BS$dates[3],
cross.midLength = lfq_results$BS$midLengths[5],
contour = TRUE,
add.values = FALSE,
agemax = 5,
hide.progressbar = TRUE)## Optimisation procuedure of ELEFAN is running.
## This will take some time.
## The process bar will inform you about the process of the calculations.# run ELEFAN with simulated annealing
res_SAbs <- ELEFAN_SA(lfq_results$BS,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 4.5,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 4,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.019744 maxTime = 30
## Emini is: -0.3184663901
## xmini are:
## 5.310221989 0.8441457965 0.2786562045 0.3598986657 0.9663543715
## Totally it used 30.019766 secs
## No. of function call is: 1190# run ELEFAN with genetic algorithm
res_GAbs <- ELEFAN_GA(lfq_results$BS,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 4.5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
# show results
RSAbs <- unlist(res_RSAbs$par)
GAbs <- unlist(res_GAbs$par)
SAbs <- unlist(res_SAbs$par)
parbs <- round(rbind(RSAbs, GAbs, SAbs),3)
parbs %>%
kbl(booktabs = T,
position="ht!",
caption = "Parametres LH to BRANSFIELD STRAIT") %>%
kable_paper("hover",
full_width = F)%>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE,
font_size=9)| Linf | K | t_anchor | C | ts | phiL | |
|---|---|---|---|---|---|---|
| RSAbs | 5.448 | 0.640 | 0.790 | 0.500 | 0.500 | 1.279 |
| GAbs | 5.297 | 0.879 | 0.428 | 0.654 | 0.516 | 1.392 |
| SAbs | 5.310 | 0.844 | 0.279 | 0.360 | 0.966 | 1.377 |
Parametrs to GERLASHE
# Response surface analyss
res_RSAgs <- ELEFAN(lfq = lfq_results$GERLASHE,
MA = 5,
Linf_range = seq(4.5, 7,
length.out = 30),
K_range = exp(seq(log(0.1),
log(2),
length.out = 30)),
C = 0.5,
ts = 0.5,
method = "cross",
cross.date = lfq_results$GERLASHE$dates[3],
cross.midLength = lfq_results$GERLASHE$midLengths[5],
contour = TRUE,
add.values = FALSE,
agemax = 5,
hide.progressbar = TRUE)## Optimisation procuedure of ELEFAN is running.
## This will take some time.
## The process bar will inform you about the process of the calculations.# run ELEFAN with simulated annealing
res_SAgs <- ELEFAN_SA(lfq_results$GERLASHE,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 4.5,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 4,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.003085 maxTime = 30
## Emini is: -0.290322014
## xmini are:
## 5.624068536 0.9826401418 0.7312699817 0.3292821422 0.6789067537
## Totally it used 30.003104 secs
## No. of function call is: 1593# run ELEFAN with genetic algorithm
res_GAgs <- ELEFAN_GA(lfq_results$GERLASHE,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 4.5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 6.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
# show results
RSAgs <- unlist(res_RSAgs$par)
GAgs <- unlist(res_GAgs$par)
SAgs <- unlist(res_SAgs$par)
pargs <- round(rbind(RSAgs, GAgs, SAgs),3)
pargs %>%
kbl(booktabs = T,
position="ht!",
caption = "Parametres LH to GERLASHE") %>%
kable_paper("hover",
full_width = F)%>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE,
font_size=9)| Linf | K | t_anchor | C | ts | phiL | |
|---|---|---|---|---|---|---|
| RSAgs | 5.448 | 0.640 | 0.790 | 0.500 | 0.500 | 1.279 |
| GAgs | 5.405 | 0.746 | 0.470 | 0.575 | 0.499 | 1.338 |
| SAgs | 5.624 | 0.983 | 0.731 | 0.329 | 0.679 | 1.492 |
Method NO LINEAR MODEL
We use firts a descomposition methot with library mixR
and then we use nlmer (Bates et al.
(2015)) to calcularte parametrs
# Agregar la nueva columna con los intervalos de longitud total
sf6 <- sf5 |>
drop_na(length_total_cm) |>
filter(ID %in% "GERLASHE",
Year==2018,
Month==5)
sf5a <- sf5 |>
drop_na(length_total_cm ) |>
filter(ID %in% "BS",
Year==2018)
# compruebo el n de comp
s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')
s_normal = select(sf5a$length_total_cm, ncomp = 2:6)
# fit a Normal mixture model
mod1 = mixfit(sf5a$length_total_cm, ncomp = 4,
pi = c(0.5, 0.5, 0.5, 0.5),
mu = c( 3, 4, 5, 6),
sd = c(0.2, 0.2, 0.2, 0.2))
mod2_weibull = mixfit(sf6$length_total_cm, family = 'weibull', ncomp = 3)
s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')
s_normal = select(sf5a$length_total_cm, ncomp = 2:6)
plot(s_weibull)
plot(s_normal)
# plot the fitted model# plot the fitted model# plot the fitted model
plot(mod1)
plot(mod2_weibull)
# fit a Normal mixture model (equal variance)
mod1_ev = mixfit(sf6$length_total_cm, ncomp = 2, ev = TRUE)References
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting
linear mixed-effects models using lme4.
Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Methot, R. D., & Wetzel, C. R. (2013). Stock
synthesis: A biological and statistical framework for fish stock
assessment and fishery management. Fisheries Research,
142, 86–99. https://doi.org/10.1016/j.fishres.2012.10.012